Advanced Programming Assignment 3: Mandelbrot zoom

Import libraries:

In [0]:
import numpy as np
import cv2 # open computer vision library, for video and image processing
import matplotlib.pyplot as plt

1. Python: mandelbrot

Define function mandelbrot with Python.

The numbers in the Mandelbrot set are all complex numbers $c$, such that the series $z_0=(0, 0)$ and $z_{t+1} = z_t^2 + c$ is bounded (i.e. it does not go to infinity). $c$ is a complex constant depending on the position in the image being processed. Set the pixel colour in relation to the number of iterations in the series before it diverges, defined as $|z_t| > 2$.

In [0]:
def mandelbrot(centre, view_side, n_max_iter = 500, view_side_pixels = 500):
  # The output is a square image, view_side is the size of the side of this
  # image in the complex plane. The view_side_pixels is the number of pixels
  # this correspond to.
  step = view_side/view_side_pixels
  # Store the pixels as a list of lists
  image = list()
  for i in range(view_side_pixels):
    image.append([])
    for j in range(view_side_pixels):
      # Find c for this pixel
      c = centre + view_side*complex(-.5, .5) + complex(j*step, -i*step)
      # Initial z
      z = complex(0, 0)
      # Iterate to generate the series
      n = 0
      # Loop until a maximum number of allowed iterations or until the 
      # series diverges.
      while n < n_max_iter and abs(z) < 2: 
        # The update
        z = z*z + c
        # Increase the iteration counter
        n += 1
      # If the series did diverge, store the number of iterations it took (our
      # colouring will be proportional to this). If not, store a zero.
      if n == n_max_iter:
        image[i].append(0)
      else:
        image[i].append(n)
  # Transform the list of lists into a numpy matrix
  return np.asarray(image)

2. Cython: mandelbrot_cy

Implement mandelbrot_cy with Cython.

load the cython extension:

In [0]:
%load_ext cython
In [5]:
%%cython -a
import numpy as np
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def mandelbrot_cy(complex centre, float view_side, int n_max_iter = 255, int view_side_pixels = 500):
  step = view_side/view_side_pixels
  cdef float step_memview = step

  cdef int i, j, n
  cdef double complex c
  cdef double complex z 
  cdef float cons_x = centre.real - view_side*.5
  cdef float cons_y = centre.imag + view_side*.5
  #cdef complex cons = centre + view_side*complex(-.5, .5)

  ret = np.zeros((view_side_pixels, view_side_pixels))
  cdef double[:, :] ret_memview = ret

  for i in range(view_side_pixels):
    for j in range(view_side_pixels):
        c = cons_x + j * step_memview + (cons_y - i * step_memview)*1j
        z = 0 + 0*1j
        n = 0
        #for n in range(n_max_iter):
         #   z = z*z + c
          #  if z.real**2 + z.imag**2 > 4:
           #     ret_memview[i, j] = n  
        while n < n_max_iter and z.real**2 + z.imag**2 <= 4: 
            z = z*z + c
            n += 1
        if n == n_max_iter:
            ret_memview[i, j] = 0
        else:
            ret_memview[i, j] = n    
  return ret
Out[5]:
Cython: _cython_magic_295e0d013b998d73940e6698bcdfbd69.pyx

Generated by Cython 0.29.15

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+01: from cython.parallel import prange
  __pyx_t_1 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+02: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 03: cimport cython
 04: 
 05: @cython.boundscheck(False)
 06: @cython.wraparound(False)
 07: @cython.cdivision(True)
+08: def mandelbrot_cy(complex centre, float view_side, int n_max_iter = 255, int view_side_pixels = 500):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_1mandelbrot_cy(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_1mandelbrot_cy = {"mandelbrot_cy", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_1mandelbrot_cy, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_1mandelbrot_cy(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  __pyx_t_double_complex __pyx_v_centre;
  float __pyx_v_view_side;
  int __pyx_v_n_max_iter;
  int __pyx_v_view_side_pixels;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("mandelbrot_cy (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_centre,&__pyx_n_s_view_side,&__pyx_n_s_n_max_iter,&__pyx_n_s_view_side_pixels,0};
    PyObject* values[4] = {0,0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        CYTHON_FALLTHROUGH;
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_centre)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_view_side)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("mandelbrot_cy", 0, 2, 4, 1); __PYX_ERR(0, 8, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (kw_args > 0) {
          PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_n_max_iter);
          if (value) { values[2] = value; kw_args--; }
        }
        CYTHON_FALLTHROUGH;
        case  3:
        if (kw_args > 0) {
          PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_view_side_pixels);
          if (value) { values[3] = value; kw_args--; }
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "mandelbrot_cy") < 0)) __PYX_ERR(0, 8, __pyx_L3_error)
      }
    } else {
      switch (PyTuple_GET_SIZE(__pyx_args)) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        CYTHON_FALLTHROUGH;
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_centre = __Pyx_PyComplex_As___pyx_t_double_complex(values[0]); if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
    __pyx_v_view_side = __pyx_PyFloat_AsFloat(values[1]); if (unlikely((__pyx_v_view_side == (float)-1) && PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
    if (values[2]) {
      __pyx_v_n_max_iter = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_n_max_iter == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
    } else {
      __pyx_v_n_max_iter = ((int)0xFF);
    }
    if (values[3]) {
      __pyx_v_view_side_pixels = __Pyx_PyInt_As_int(values[3]); if (unlikely((__pyx_v_view_side_pixels == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
    } else {
      __pyx_v_view_side_pixels = ((int)0x1F4);
    }
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("mandelbrot_cy", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 8, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_295e0d013b998d73940e6698bcdfbd69.mandelbrot_cy", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_mandelbrot_cy(__pyx_self, __pyx_v_centre, __pyx_v_view_side, __pyx_v_n_max_iter, __pyx_v_view_side_pixels);

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_mandelbrot_cy(CYTHON_UNUSED PyObject *__pyx_self, __pyx_t_double_complex __pyx_v_centre, float __pyx_v_view_side, int __pyx_v_n_max_iter, int __pyx_v_view_side_pixels) {
  PyObject *__pyx_v_step = NULL;
  float __pyx_v_step_memview;
  int __pyx_v_i;
  int __pyx_v_j;
  int __pyx_v_n;
  __pyx_t_double_complex __pyx_v_c;
  __pyx_t_double_complex __pyx_v_z;
  float __pyx_v_cons_x;
  float __pyx_v_cons_y;
  PyObject *__pyx_v_ret = NULL;
  __Pyx_memviewslice __pyx_v_ret_memview = { 0, 0, { 0 }, { 0 }, { 0 } };
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("mandelbrot_cy", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_6);
  __PYX_XDEC_MEMVIEW(&__pyx_t_7, 1);
  __Pyx_AddTraceback("_cython_magic_295e0d013b998d73940e6698bcdfbd69.mandelbrot_cy", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XDECREF(__pyx_v_step);
  __Pyx_XDECREF(__pyx_v_ret);
  __PYX_XDEC_MEMVIEW(&__pyx_v_ret_memview, 1);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__19 = PyTuple_Pack(15, __pyx_n_s_centre, __pyx_n_s_view_side, __pyx_n_s_n_max_iter, __pyx_n_s_view_side_pixels, __pyx_n_s_step, __pyx_n_s_step_memview, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n, __pyx_n_s_c, __pyx_n_s_z, __pyx_n_s_cons_x, __pyx_n_s_cons_y, __pyx_n_s_ret, __pyx_n_s_ret_memview); if (unlikely(!__pyx_tuple__19)) __PYX_ERR(0, 8, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__19);
  __Pyx_GIVEREF(__pyx_tuple__19);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_295e0d013b998d73940e6698bcdfbd69_1mandelbrot_cy, NULL, __pyx_n_s_cython_magic_295e0d013b998d7394); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 8, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_mandelbrot_cy, __pyx_t_1) < 0) __PYX_ERR(0, 8, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_codeobj__20 = (PyObject*)__Pyx_PyCode_New(4, 0, 15, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__19, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_root_cache_ipython_cython__cyth, __pyx_n_s_mandelbrot_cy, 8, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__20)) __PYX_ERR(0, 8, __pyx_L1_error)
+09:   step = view_side/view_side_pixels
  __pyx_t_1 = PyFloat_FromDouble((__pyx_v_view_side / ((float)__pyx_v_view_side_pixels))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_step = __pyx_t_1;
  __pyx_t_1 = 0;
+10:   cdef float step_memview = step
  __pyx_t_2 = __pyx_PyFloat_AsFloat(__pyx_v_step); if (unlikely((__pyx_t_2 == (float)-1) && PyErr_Occurred())) __PYX_ERR(0, 10, __pyx_L1_error)
  __pyx_v_step_memview = __pyx_t_2;
 11: 
 12:   cdef int i, j, n
 13:   cdef double complex c
 14:   cdef double complex z
+15:   cdef float cons_x = centre.real - view_side*.5
  __pyx_v_cons_x = (__Pyx_CREAL(__pyx_v_centre) - (__pyx_v_view_side * .5));
+16:   cdef float cons_y = centre.imag + view_side*.5
  __pyx_v_cons_y = (__Pyx_CIMAG(__pyx_v_centre) + (__pyx_v_view_side * .5));
 17:   #cdef complex cons = centre + view_side*complex(-.5, .5)
 18: 
+19:   ret = np.zeros((view_side_pixels, view_side_pixels))
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_view_side_pixels); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_view_side_pixels); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_6 = PyTuple_New(2); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_GIVEREF(__pyx_t_3);
  PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_3);
  __Pyx_GIVEREF(__pyx_t_5);
  PyTuple_SET_ITEM(__pyx_t_6, 1, __pyx_t_5);
  __pyx_t_3 = 0;
  __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  __pyx_t_1 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_5, __pyx_t_6) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_6);
  __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_ret = __pyx_t_1;
  __pyx_t_1 = 0;
+20:   cdef double[:, :] ret_memview = ret
  __pyx_t_7 = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(__pyx_v_ret, PyBUF_WRITABLE); if (unlikely(!__pyx_t_7.memview)) __PYX_ERR(0, 20, __pyx_L1_error)
  __pyx_v_ret_memview = __pyx_t_7;
  __pyx_t_7.memview = NULL;
  __pyx_t_7.data = NULL;
 21: 
+22:   for i in range(view_side_pixels):
  __pyx_t_8 = __pyx_v_view_side_pixels;
  __pyx_t_9 = __pyx_t_8;
  for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
    __pyx_v_i = __pyx_t_10;
+23:     for j in range(view_side_pixels):
    __pyx_t_11 = __pyx_v_view_side_pixels;
    __pyx_t_12 = __pyx_t_11;
    for (__pyx_t_13 = 0; __pyx_t_13 < __pyx_t_12; __pyx_t_13+=1) {
      __pyx_v_j = __pyx_t_13;
+24:         c = cons_x + j * step_memview + (cons_y - i * step_memview)*1j
      __pyx_v_c = __Pyx_c_sum_double(__pyx_t_double_complex_from_parts((__pyx_v_cons_x + (__pyx_v_j * __pyx_v_step_memview)), 0), __Pyx_c_prod_double(__pyx_t_double_complex_from_parts((__pyx_v_cons_y - (__pyx_v_i * __pyx_v_step_memview)), 0), __pyx_t_double_complex_from_parts(0, 1.0)));
+25:         z = 0 + 0*1j
      __pyx_v_z = __Pyx_c_sum_double(__pyx_t_double_complex_from_parts(0, 0), __Pyx_c_prod_double(__pyx_t_double_complex_from_parts(0, 0), __pyx_t_double_complex_from_parts(0, 1.0)));
+26:         n = 0
      __pyx_v_n = 0;
 27:         #for n in range(n_max_iter):
 28:          #   z = z*z + c
 29:           #  if z.real**2 + z.imag**2 > 4:
 30:            #     ret_memview[i, j] = n  
+31:         while n < n_max_iter and z.real**2 + z.imag**2 <= 4:
      while (1) {
        __pyx_t_15 = ((__pyx_v_n < __pyx_v_n_max_iter) != 0);
        if (__pyx_t_15) {
        } else {
          __pyx_t_14 = __pyx_t_15;
          goto __pyx_L9_bool_binop_done;
        }
        __pyx_t_15 = (((pow(__Pyx_CREAL(__pyx_v_z), 2.0) + pow(__Pyx_CIMAG(__pyx_v_z), 2.0)) <= 4.0) != 0);
        __pyx_t_14 = __pyx_t_15;
        __pyx_L9_bool_binop_done:;
        if (!__pyx_t_14) break;
+32:             z = z*z + c
        __pyx_v_z = __Pyx_c_sum_double(__Pyx_c_prod_double(__pyx_v_z, __pyx_v_z), __pyx_v_c);
+33:             n += 1
        __pyx_v_n = (__pyx_v_n + 1);
      }
+34:         if n == n_max_iter:
      __pyx_t_14 = ((__pyx_v_n == __pyx_v_n_max_iter) != 0);
      if (__pyx_t_14) {
/* … */
        goto __pyx_L11;
      }
+35:             ret_memview[i, j] = 0
        __pyx_t_16 = __pyx_v_i;
        __pyx_t_17 = __pyx_v_j;
        *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_ret_memview.data + __pyx_t_16 * __pyx_v_ret_memview.strides[0]) ) + __pyx_t_17 * __pyx_v_ret_memview.strides[1]) )) = 0.0;
 36:         else:
+37:             ret_memview[i, j] = n
      /*else*/ {
        __pyx_t_18 = __pyx_v_i;
        __pyx_t_19 = __pyx_v_j;
        *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_ret_memview.data + __pyx_t_18 * __pyx_v_ret_memview.strides[0]) ) + __pyx_t_19 * __pyx_v_ret_memview.strides[1]) )) = __pyx_v_n;
      }
      __pyx_L11:;
    }
  }
+38:   return ret
  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(__pyx_v_ret);
  __pyx_r = __pyx_v_ret;
  goto __pyx_L0;

3. Compare Execution Times

In [0]:
import timeit
In [7]:
%timeit mandelbrot(complex(-.74303, .126433), .01611, view_side_pixels = 200)
%timeit mandelbrot_cy(complex(-.74303, .126433), .01611, view_side_pixels = 200)
1 loop, best of 3: 1.76 s per loop
10 loops, best of 3: 28.5 ms per loop

4. Generate Image and Video

4.1. Image

View the generated image.

A colourmap is applied to the raw counts of iterations before diverging. Note that opencv use BGR instead of RGB colour ordering.

In [8]:
# Make the image
image = mandelbrot_cy(complex(-.5, 0), 2, view_side_pixels = 200)
# image = mandelbrot(complex(-.74303, .126433), .01611, view_side_pixels = 200)

# Apply colourmap
img = np.asarray(255*image/np.max(image), dtype=np.uint8)
img = cv2.applyColorMap(img, cv2.COLORMAP_JET)

# Optional resize (makes it look better at higher resolutions)
# img = cv2.resize(img, (img.shape[0]//2, img.shape[1]//2), interpolation=cv2.INTER_CUBIC)

# Write the image and view it
cv2.imwrite('fractal.png', img)
from IPython.display import Image
Image('fractal.png')
Out[8]:

4.2. Video

Generate a 500-frame video at a 1000x1000 pixel resolution. However, since the video is too large, we only show every 25 frame:

In [0]:
def framemaker(zoom):
  # Create the image for the given parameters
  image = mandelbrot_cy(complex(-0.743643887037158704752191506114774, 
                                 0.131825904205311970493132056385139), view_side=zoom, n_max_iter=255, view_side_pixels=1000)
  # Apply colourmap
  img = np.asarray(255*image/np.max(image), dtype=np.uint8)
  img = cv2.applyColorMap(img, cv2.COLORMAP_JET)
  # Optional resize (makes it look better at higher resolutions)
  img = cv2.resize(img, (img.shape[0]//2, img.shape[1]//2), interpolation=cv2.INTER_CUBIC)
  return img
In [0]:
from multiprocessing import Pool
import os

def parallizer():
  zoom = np.exp(np.linspace(np.log(1), np.log(.00001), 500))
  with Pool(processes=os.cpu_count()) as pool:
    frames = pool.map(framemaker, zoom)
  return frames
In [0]:
frames = parallizer()
In [12]:
# Show every 25 rendered frame
from IPython.display import Image, display
Image('fractal.png')
for i in range(0,500,25):
  cv2.imwrite('fractal.png', frames[i])
  display(Image('fractal.png'))

Comments

  1. In the function mandelbrot_cy, when updating the matrix ret (of shape (view_side_pixels, view_side_pixels)), because there are two for loops, one nested within another, which makes it difficult to implement map to update the matrix in parallel. In other words, because the update process of the matrix needs to remember where it was in the previous step (i.e. position denoted by (i ,j)), therefore it cannot be run in parallel.

  2. We tried to use prange as an alternative to map to let the outermose for-loop be run in parallel. However, this slowed down our function by maybe a quarter. This is in part because we needed to use a for-loop as the functions innermost loop (see the commented out code in mandelbrot_cy) which was slower than our finally used while-loop.

  3. After speeding up, our implementation with Cython is around 60 times than the original implementation using Python.